#FH_maps.R
#Replication code for maps in 
#Feigenbaum and Hall
#How Legislators Respond To Localized Economic Shocks: Evidence From Chinese Import Competition

#load libraries
library(maptools)
library(ggplot2)
library(foreign)
library(data.table)
library(RColorBrewer)
library(grid)

#load shapefiles from Lewis et al (111th and 106th)
cd106 <- readShapePoly("districts106.shp")
cd111 <- readShapePoly("districts111.shp")

#load cz shapefile
cz <- readShapePoly("CZ.shp")

#prep for ggplot2
cd106.df <- fortify(cd106)
cd106.df <- data.frame(cd106.df, cd106@data[cd106.df$id, ])

#prep for ggplot2
cd111.df <- fortify(cd111)
cd111.df <- data.frame(cd111.df, cd111@data[cd111.df$id, ])

#prep for ggplot2
cz.df <- fortify(cz)
cz.df <- data.frame(cz.df, cz@data[cz.df$id, ])

#trade shock data
tradeshock <- as.data.table(read.dta("fh_final_analysis.dta"))
tradeshock <- tradeshock[, list(state, dist, decade, x, z)]

#convert state abbrevs to state names
tradeshock$STATENAME <- state.name[match(tradeshock$state, state.abb)]

#DISTRICT is a factor, needs to be numeric
cd106.df$dist <- as.numeric(as.character(cd106.df$DISTRICT))
cd111.df$dist <- as.numeric(as.character(cd111.df$DISTRICT))

#dist = 0 for at large districts in the map data
cd106.df$dist <- cd106.df$dist * (cd106.df$dist != 0) + 1 * (cd106.df$dist == 0)
cd111.df$dist <- cd111.df$dist * (cd111.df$dist != 0) + 1 * (cd111.df$dist == 0)

# Restrict to continental USA
cd106.df <- cd106.df[cd106.df$long < -50 & cd106.df$long > -135 & cd106.df$lat < 50 & cd106.df$lat > 22, ]
cd111.df <- cd111.df[cd111.df$long < -50 & cd111.df$long > -135 & cd111.df$lat < 50 & cd111.df$lat > 22, ]
cz.df <- cz.df[cz.df$long < -50 & cz.df$long > -135 & cz.df$lat < 50 & cz.df$lat > 22, ]

#tradeshocks for decade 1 merge with 106th CD, decade 2 merge with 111th CD
ts.106 <- tradeshock[decade == 1,]
ts.111 <- tradeshock[decade == 2,]

#Figure 3
#CZ and CDs in 2000s

pdf("Figure_3.pdf",width=7,height=5)

mapPlot <- ggplot() + theme_bw() +
  theme(
    plot.background = element_blank()
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.border = element_blank()
    ,axis.title = element_blank()
    ,axis.text = element_blank()
    ,axis.ticks = element_blank()
    ,legend.position="bottom"
    ,legend.key.width = unit(2, "cm")
  )

mapPlot <- mapPlot + geom_polygon(data = cz.df, aes(x = long,
                                  y = lat,
                                  group = group),
                                  fill = NA,
                                  colour = "gray85",
                                  lwd = 1/5)

mapPlot <- mapPlot + geom_polygon(data = cd111.df, aes(x = long,
                                  y = lat,
                                  group = group),
                                  fill = NA,
                                  colour = "black",
                                  lwd = 1/5)


# Make a conical projection
# -- otherwise it looks funny
mapPlot <- mapPlot + coord_map(project="conic", lat0 = 30)
plot(mapPlot)
dev.off()

rm(mapPlot)

#Figure 4
#trade shocks is 2000-07

cd111.df <- merge(cd111.df, ts.111, by = intersect(names(cd111.df), names(ts.111)))
cd111.df <- cd111.df[with(cd111.df, order(STATENAME,dist,order)),]

myPalette <- colorRampPalette(brewer.pal(7, "Reds"))

pdf("Figure_4.pdf",width=7,height=5)

mapPlot <- ggplot(cd111.df) + theme_bw() +
  theme(
    plot.background = element_blank()
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.border = element_blank()
    ,axis.title = element_blank()
    ,axis.text = element_blank()
    ,axis.ticks = element_blank()
    ,legend.position="bottom"
    ,legend.key.width = unit(1, "cm")
    ,legend.text = element_text(size = 8)
    )

mapPlot <- mapPlot + geom_polygon(aes(x = long,
               y = lat,
               group = group,
               fill = x),
               colour = "white",
               lwd = 1/100)
# Make a conical projection
mapPlot <- mapPlot + coord_map(project="conic", lat0 = 30)
mapPlot <- mapPlot + scale_fill_gradientn("Trade Shock", colours = myPalette(100), trans = "sqrt")
plot(mapPlot)
dev.off()

rm(mapPlot)

#Figure 8
#CZs and CDs in 1990s (a) and 2000s (b)

pdf("Figure_8_a.pdf",width=7,height=5)

mapPlot <- ggplot() + theme_bw() +
  theme(
    plot.background = element_blank()
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.border = element_blank()
    ,axis.title = element_blank()
    ,axis.text = element_blank()
    ,axis.ticks = element_blank()
    ,legend.position="bottom"
    ,legend.key.width = unit(2, "cm")
  )

mapPlot <- mapPlot + geom_polygon(data = cz.df, aes(x = long,
                                  y = lat,
                                  group = group),
                                  fill = NA,
                                  colour = "gray85",
                                  lwd = 1/5)

mapPlot <- mapPlot + geom_polygon(data = cd106.df, aes(x = long,
                                  y = lat,
                                  group = group),
                                  fill = NA,
                                  colour = "black",
                                  lwd = 1/5)


# Make a conical projection
# -- otherwise it looks funny
mapPlot <- mapPlot + coord_map(project="conic", lat0 = 30)
plot(mapPlot)
dev.off()

rm(mapPlot)

pdf("Figure_8_b.pdf",width=7,height=5)

mapPlot <- ggplot() + theme_bw() +
  theme(
    plot.background = element_blank()
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,panel.border = element_blank()
    ,axis.title = element_blank()
    ,axis.text = element_blank()
    ,axis.ticks = element_blank()
    ,legend.position="bottom"
    ,legend.key.width = unit(2, "cm")
  )

mapPlot <- mapPlot + geom_polygon(data = cz.df, aes(x = long,
                                  y = lat,
                                  group = group),
                                  fill = NA,
                                  colour = "gray85",
                                  lwd = 1/5)

mapPlot <- mapPlot + geom_polygon(data = cd111.df, aes(x = long,
                                  y = lat,
                                  group = group),
                                  fill = NA,
                                  colour = "black",
                                  lwd = 1/5)


# Make a conical projection
# -- otherwise it looks funny
mapPlot <- mapPlot + coord_map(project="conic", lat0 = 30)
plot(mapPlot)
dev.off()

rm(mapPlot)
